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CALCULATION OF EPHEMERIDES FROM INITIAL VALUES 


by 

Karl Stumpff 

Goddard Space Flight Center 


SUMMARY 

A short report is presented on a method for calculating un- 
disturbed ephemerides (coordinates of location and velocity) of a 
planet or a comet when its initial time values are given. Thus far 
this method has been available only in German language publica- 
tions. The orbital elements need not be known for this method, 
which is applicable without formal variations for all types of or- 
bits. In particular, the singularity of classical methods is avoided 
by the transfer from elliptical to hyperbolic orbits. In place of 
Kepler 1 s equation, a transcendent main equation appears which is 
valid for all types of orbits and becomes rational for circular 
and parabolic orbits. The formulas for the calculation of location 
and velocity coordinates are simple and especially well suited 
for electronic computers. The optimal area for application (small 
and medium intermediate times) coincides with the requirements 
for orbit determination, orbit correction, and the calculation of 
special perturbations. 
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CALCULATION OF EPHEMERIDES FROM INITIAL VALUES 


by 

Karl Stumpff* 

Goddard Space Flight Center 


THE PROBLEM 

The undisturbed orbit of a celestial body orbiting around the sun is definitely deter- 
mined if the coordinates of its location and velocity are given in relation to the sun. From 
these six quantities, x 0 , y 0 , z 0 , x Q , y Q , z Q — the local elements * of the orbit — the six clas- 
sical orbital elements, i , fi, o>, a, e , T, can be derived. And then x, • • , z can be computed 
for any other given time t with the help of the orbital elements. This process requires the 
solution of the transcendental Kepler equation for every t if the orbit is an ellipse with a 
relatively small eccentricity. For parabolas, near-parabolic ellipses and hyperbolas, as 
well as for greatly eccentric hyperbolas (e » i), other relationships occur in place of the 
Kepler equation, so that very different computational processes must be used for these 
orbits. 

The possibility of arriving at a solution without using the classical orbit elements and 
without variations owing to individual case differences is indicated by the Taylor develop- 
ment of coordinates; for instance 


x( T ) 


+ X q T 


+ 2 x o T 


x‘ r 3 + 


(1) 


where r = k (t -tj denotes the "intermediate time" expressed in units of 1/k = 58 < ?13244, 
and the higher derivatives x' 0 , x‘ 0 , • • • , valid for t Q , can be derived from the equations of 
motion 


*+-7 = 0, •••*+-7 = 0 (r* = X* + y2 + **) . (2) 

r* 5 r J 

But these formulas which are independent of orbit shape are only sufficiently convergent 
for small intermediate times r, and thus are utilized only occasionally for first orbit 


♦NAS-NASA Research Associate; Professor Emeritus, Gottingen University. 

tThe local elements of the orbit are so named because their values depend upon the location of the orbit. 
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determination. The method described herein shows that it is also possible to develop ex- 
act closed formulas where the coordinates x(r) • • • z(r) can be represented directly as 
functions of the local elements and an intermediate time of any given length, r . 

THE LOCAL INVARIABLES 

Two-body motion occurs in a plane given by the two 
vectors p 0 (x 0 , y 0 , z 0 ) , p 0 (x 0 , y 0 , z 0 ) , which define the lo- 
cation and velocity of the planet at the time t 0 . With the 
exception of the situation wherein p 0 and p 0 are coinci- 
dental or oriented opposite to one another (and this case, 
which indicates straight-line motion either directly toward 
or away from the sun, shall be discounted). Since p(r) , 
and also p(r) , can be broken down into components ac- 
cording to the directions p 0 and p 0 , we have 

p(r) = p 0 F + p 0 G, P (t) = p 0 F + p 0 G , (3) 

where F and G are functions of the local elements and the intermediate time r. 

Further, it is clear that the relative locations of the vector P (r) and the scalar 
quantities F and G for any given value of r, depend only upon the geometric character- 
istics of Figure 1, defined by p 0 and p 0 — in other words, by the values r 0 or V Q of these 
vectors and the angle 5 between them — and not upon the orientation of this figure in an 
area, or space, coordinate system. From this, however, it follows that F and G contain 
the local elements only in such form as 

x o 2 + y 0 2 + 2 o 2 = (P0P0) = p n = r o 2 . 

x o x o + y 0 y 0 + z o i o = (P 0 P 0 ) = p i 2 = r o v o cosS 

x 0 2 + + *<? = (P0P0) = P22 = v 0 2 ■ 

which we call local invariables, because they are independent of the coordinate system 
selected, although they vary from point to point of the orbit. In actuality, if we set 


then, according to Equation 2, 




Figure 1 — Geometric repre- 
sentation of two-body motion 
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X - “ X/X , ^ 

X = -x£ - x/i , > (6) 

x iv = -xjll - 2x/i “x/x = x(/x 2 - /x) - 2xA . etc.; > 

and we find that, in the development of Equation 1, all higher derivations of x can be rep- 
resented linearly through x and x, with coefficients which are dependent only upon /x,/x, 


Again, these quantities can be expressed through the invariables (Equation 4) as can 
be demonstrated easily through step by step differentiation of Equation 5. We introduce 
the quantities 





Cr 


r 

r 


xx + yy + zz 
r 2 



r (?) 




x 2 + y 2 + z 2 _ p 22 

~ 2 “ p^T 


as the fundamental invariables . By differentiating them and substituting -fix, • • • -/xz for 
x, • • • z , we find the relations: 


/t - - 3/x o* , or f 


a - a; - /x - 2 cr 2 ; 


^ ( 8 ) 


a; = — 2cr{ fj. + co) ; 


and also 


ju - 15/xcr 2 “3/xoj + 3/x 2 ; 

Ju* = - 105 /x cr 3 + 45 ficr co - 42 /x 2 cr ; 


(9) 


etc. 


Thus, all derivations of /x (or r), cr, ^ are themselves functions of these quantities. It fol- 
lows that F and G are functions of the intermediate time r and of /x 0 , cr Q , w 0 , if these are the 
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values of the invariables for t = t 0 (r = o). Finally, from Equation 1 (when Equations 6 
and 9 are introduced, we see that p and cj always occur in connection with r 2 , and that a al- 
ways occurs in connection with r, and that F and (l/r)G can be written as power series of 
, cr 0 T, and oj 0 t 2 . 

However, experience has shown that it is advantageous to introduce, in addition to 
Equation 7, other complex invariables such as 

e = co - p, p - 2p - co, & = co - cr 2 (10) 

which have a special meaning in two-body motion, and can replace one or the other of the 
invariables (Equation 7) if necessary. The derivatives of Equation 10 are 

k = - cr ( 2co ~~ p) , p = “ 2pcr, & - - 4tfcr . (11) 

In particular, by the elimination of cr from r = rcr, p - - 2 pa, $ = -4tfa, we obtain 



whose integration yield the integrals of energy and areas: 

r 2 p = constant = , and r 4 t? = constant = p = a(l - e 2 ) . (12) 

It follows from this that the quantity p serves as a criterion for the orbit shape, since 
P> ofor ellipses (a>0), p = 0 for parabolas (a = ®) and p< 0 for hyperbolas (a<o). 

All geometric quantities of a conic-section orbit can be expressed in more or less 
simple form by means of the invariables (Equations 7 and 10); for example, we readily 
find that 


e cos v 


p - _ t? i „ 

y “ 1 - — - 1 ; e cos E 


1 a p ’ 


. cr r s in v _ o' Vp 

e sin v - r fp “ 5 e sinE " — ^==~~ “ — p — 


V* ; 


r as) 


E = ^ ; 


i2 = 


€ 2 + per 2 


■! = ^ 
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Finally, we obtain 


r = r a , 



from which, by the elimination of 0 - and e to- 2 , there follows 


r* + 



0 . 




[ (14) 


(15) 


INTRODUCTION OF A NEW INDEPENDENT VARIABLE 


In order to integrate Equation 15, we introduce a new variable q in place of the time r 


dr - rdq 


f 

•'n 


r ( q ) dq . 


For the initial epoch r = o, and also q = 0; since 1/r is essentially positive as long as r 
remains finite, q increases proportionally in the same sense as r, so that q(r) represents 
a definitely reversible function. Therefore, if we set* 



then Equation 15 becomes 


(17) 


Derivations with respect to q are indicated by primes. 
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However, 



thus a 2 = (l -r")/r is constant and Equation 17 assumes the simple form 

r"'+a 2 r' = 0 ( 18 ) 

whose complete integral, with three constants a, b, c is, in closed form, 


a+bcosaq + c sinaq 

(for 

a 2 > 0) , 



a + bq + cq 2 

(for 

a 2 = 0) , 

> 

(19) 

a + b cosh /Bq + c sinh/ 3 q 

(for 

~/ 3 2 = a 2 < 0) , 

7 J 




The integrals (Equation 19) still contain the case variations of the classic two-body theory. 
Our requirements for formulas applicable without difference to all orbit shapes can be 
satisfied by utilizing the following expedient: 

The differential Equation 18 has the particular integral 


r - cos aq . 

By continued integration with respect to q, a series of auxiliary functions are obtained: 

c n = f f cosaf(d<f) n . (20) 

q Jo Jo 

The beginning of this series of auxiliary functions (which we shall term c-functions) is: 

sin aq 1 - cos aq _ aq - sin aq 

c 0 = cosaq; c, = —55- ! c, = ; c 3 — (21) 

By setting ^ = aq, we can also write the c-functions in the form of the constantly con- 
vergent power series 


c> 2 ) 


1 k 2 \ 4 \ 6 

nT~ ( n + 2 ) ! + ( n + 4 ) ! _ ( n + 6 ) ! + 


( 22 ) 
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For this there applies the differential relationship 


35 ( c «+i q " +1 ) " c « q " 


(23) 


and the recurrence formula 


7JT = C n + ^ C n + 2 


(24) 


Now we can write r(q) as a Taylor series: 


r(q) - 'o+Vq+i r 0 V+ir»q* + 


or, by substituting Equation 24 for the reciprocal factorials 


r ( 


q) = r 0 + V ( c i + k2c 3 ) q + r o" ( C 2 + ^ C 4 ) q2 + r o”( C 3 + ^ C S ) 


If we again substitute k 2 = a 2 q 2 and arrange the preceding equation according to powers of 
q, we obtain 

r(q) = r Q + c,r 0 'q + c 2 r 0 " 2 + c 3 (r 0 '» + a 2 r 0 ') q 3 + c 4 ( t ™ + a 2 r 0 ") q 4 + ••• 

in which all terms of third and higher orders disappear because of Equation 18. 

We now have the closed expression 


r(q) - r 0 + C^’q + c 2 r 0 "q 2 


(25) 


But now, because of Equations 12, 14, and 16, 


r 

q 


t" -4" (r 2 + r r ) = r 3 e , 


1 -r* _ 1 - r 3 e _ 1 -r 3 ( M ~/o) 


t 2 P = -=■ 
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we obtain 


k 2 = a 2 q 2 = — = p( rq } 2 

as the argument of the of unctions. The closed expression of Equation 25 
written 

r(q) = r 0 [l + Cl a 0 (r 0 q) + c 2 e 0 (r 0 q)2] , 

where c n = c n [/> 0 (r 0 q) 2 ] . 

Finally, if we set 

£ 0 = ^0 r2 > % = a o T - to = € o t2 > *o = P 0 t2 


and, in place of q, introduce the variable 

2 


r o q 


we obtain, instead of Equation 26: 

r = r 0 (l +Cj V 0 z +c 2 £ 0 z 2 ) , 


where c n = c n ( Xo z 2 ); or 

r = r o( c o +c i 7 ’o z + c 2^o z2 ) - 
where c 0 = l -c 2Xo z 2 and Xo = i 0 -C 0 . 

THE MAIN EQUATION 

To clarify the relationship between z andr, consider that 

dr = rdq = % (l + Cl q • r Q <7 0 +c 2 q 2 • r fl 2 €„) dq . 

Upon integrating and setting 

Jc n q n dq = c n + 1 q n + I 


can now be 

(26) 


(27) 


(28) 


(29a) 


(29b) 


( 30 ) 
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according to Equation 23, Equation 30 becomes 

T = r 0 C l + C 2 r 0 2 ° r 0 < 5 2 + C 3 r 0 3e 0 < 5 3 • 

If we divide this by r and substitute z (Equation 28), we finally obtain 

1 = z + c 2 cr 0 Tz 2 + c 3 e 0 r 2 z 3 

or 


1 - z + c 2 t}q z 2 + c 3 z 3 ; 


f ° r C n = C n(*0 z2 ) • 


(31a) 


If we add 0 = + c 3 x 0 z 2 -l) to the above equation, where x 0 = € 0 ~i 0 > it can also be 

written: 

1 = Cj Z + C 2 7) 0 Z 2 + C 3 ^ 0 Z 3 , (31 b ) 

This transcendental equation represents the Kepler equation as developed here and 
the analogous equations for parabolic and hyperbolic motion. The Equations 29 and 31 are 
valid for all forms of orbits: in hyperbolic motion, the argument of the c-functions becomes 
negative, since p 0 and y Q - p Q r 2 become negative, but the functions themselves remain 
real. In the case of the parabola, the c n are approaching their constant terml/n!; thus 
\ 2 = x 0 z2 becomes zero; and the main equation becomes a rational third degree equation: 

1 = z + T % z2 + T £o z3 • ( 32 ) 

For circular orbits {e = o), o- Q = e Q = o from Equation 13, and v 0 = = therefore 

Equation 31 has the trivial form z = 1. 

Equation 31a, termed the main equation of the two-body problem, can also be developed 
directly from the Kepler equation. If we write t Q and t for two times, 

t 

it-. 
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where T is the time of passage through the perihelion, and then set k = E -e 0 , the differ- 
ence of both equations is 


k - e [sin (E q +\) - sin eJ 



or 


(l -e cos E Q ) sin\ + e sinE Q ( 1 -cos \) + (V-sin>J = — — — • 

By substituting in the foregoing the expressions of Equation 13 for e cosE 0 , e sinE Q3 
a" 3/2 , and the c-f unctions 

c i(^ 2 ) = sin c 2 (\ 2 ) = (1-cos X.)A 2 , and c 3 (k 2 ) = sin k)/k 3 , 

we obtain 


p 0 jfT 0 

c f ~~r . — k + c -O' — — X 2 + cA 3 


Tp 0 fo 


'* ^0 


' 2 0 /“o 


If we then divide through by term on the right side and put 


: 1 = 1 _ ^ C 3- X = ZT /V a 0 T ~ V S 0 T * ~ £o’ 


we obtain the main equation in the form of Equation 31. The quantity k = E-E 0 = zr 
introduced in the argument of the c-functions, is real only for positive p 0 or for elliptical 
orbits, and is identical with k = a q (which we used in the previous section) since 
a = 1/fr = r 0 yz^ and k = r 0 q yp^ = zr^,from Equation 28. 


From z = (r o /V) q and q = 



it follows that 


z 



or, if we set 


that 


r = r 0 A 


r o ( 1 +C 1 7 ?0 Z+ c 2 ^0 z2 ) 



dr 

A 


(33) 


Z 


T 


(34) 
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Therefore z is the mean value for l/A = r Q /r over the time interval (0, r) and is always 
a positive quantity. For circular orbits r = r 0 , and therefore A = l from Equation 33 and 
z = 1 from Equation 34, which values we could also ascertain from the main equation. 

The form of the main equation valid for parabolic orbits (Equation 32) can be derived 
from the well known cubic equation for the tangent of the half true anomaly, which replaces 
the Kepler equation for e = l: 


v 1 v t — T 

tan ~2 + -3 tan 3 ~ 2k 


where p = the parameter of the parabola. By subtracting 


v o 1 , v o 

tan ~2 + -3- tan J 


= 2k 


*0 - T 


(for t = t 0 ) from the above equation and by letting 


we obtain the relation 


tan ~2 ~ tan 


B - 1 + tan ~2 tan 


(35) 


+4 A2 ) 


2 T 


(36) 


Since, from Equation 13, 


tan v - 


2 tan T _ 


1 - tan 2 "2" M 


0 V t?-U V 

tan T + 2-^-^-tan T 


1 , 


we deduce as a solution of this quadratic equation 

: " * ± /^ 2 - * ( 2 m - * " ct2 ) ] 


tan 



mi n ii 1 1 in mi i 


ii ■■■■hi mi n 


According to Equation 10, 


2/2 - t ) - cr 2 = (2/2—0)) + co-cr 2 - d - P ; 


but in the case of the parabola p - 0, either 


v _ 2 


tan -Ty - 


yw •/& 


Here the first of these solutions applies, since o- and v disappear in the perihelion. 
Now, from Equation 16 r' = dr/dq = r 2 a ) r " = r 3 e; and from Equation 25 

r ' = r o' (vO* + r o" MT 


r o 2cr o c o‘l 0 + r o 3e o( c i ( j) 


3 r o 2 ( c o CT o + c i e o T o*) 


- r o 2 ( c o CT o + c x 6 0 ZT ) ■ 

Therefore in the case of the parabola ( c 0 “ c 1 =l,/o=/ 2 -€=oj , we have 


i / v 

— r (cr + £ a zt) 
A2 l 0 0 / 


^0 + Vq zT 
A 2 


Furthermore r 4 tf - r 0 4 tf 0 from Equation 12, thus 


Since 


2 _ yt: 
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and 


v 

tan ~ 


CT 0 + M o ZT 


& ~ /r 0 


we find from Equation 35 that 


A = 


M 0 zr 


/ \ 
B = 1 + *7 K + M 0 ZT ) 


M 0 , X 

*7 ( 2 +or o ZT ) 


Now, because p = 0, we may set 


r 0 2 + *0 = V 0 


Equation 36 then has the form 


^o z 




1 Mo 2 * 2 t2 


-^-(2+a fl z r) +3 


2r 


or 


1 9 1 _ , _ 1 

z + 2" % z + ~6^o z ~ ~~2 

r 0 




(37) 


where we let = y^r 2 , r) 0 - <x^r according to Equation 27. This equation is identical with 
Equation 32; however, 


p 

2 


1 + tan 2 ~~2 



r o*o 
*0 + a o 


r o *, 

2 M 0 


so that 




and therefore the right side of Equation 37 equals unity. Furthermore, for parabolas, 

m 0 = e 0 ; thus = c«- 
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F AND G AS FUNCTIONS OF T 

If, with the help of the main equation, z has been determined as a function of r and of 
the local invariables of the initial epoch, the quantities F and G as well as their derivatives 
F and G can be expressed easily as functions of z. If in the vectorial differential equation 


P = -MP 


we set 


P = P 0 F + Po G * P = Po F + Po G 
in the manner of Equation 3, there follows that the identity 

p 0 (F + ijF) + p 0 (G + pG) = 0 


is fulfilled only if 


F + fJF 


0 , G + fjG = 0 . 


(38) 


If a new variable q is introduced (with q = l/r as above), we then have 


F = 






r (39) 


since r = rcr; and we obtain 

rF" - r 2 oF' + F = 0 

instead of Equation 38, If we differentiate the above equation again with respect to q, we 
obtain 

rF'" + r'F" - 2rr'oF' - r 2 cr'F' - rVF" + F # = 


0 
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or, since 


cr » 

ex' - — - r<x = r (e ~ 2^ 2 ) , 

q v 

1 - r 3 /j, , 


rF'" + r 3 (/x-e)F' - r(F"' + r 2 pF') = 0 . 


But r 2 p = a 2 is constant; thus 


F"' + a 2 F f = 0 , 


and 


G /m + a 2 G f = 0 , 


which means that the same differential equation will suffice for F(q) and G{q) as for r(q) . 
Therefore, since exactly the same operations can be used on these two functions as on 
r ( q) in Equation 25, we have 


F(q) = F q + Cl F 0 ' q + c 2 F 0 "q 2 
G(q) = G q + Cl G 0 ' q + c.G/q 2 . 


(40) 


Now from 


x(t) = x 0 F(r) + x q G(t) , 
x(r) = x Q F(t) + x Q G(t } 

it follows that, for r = 0, 


F o 

= 1, 

F o 

= 0 

G o 

= o, 

G o 

= 1 



16 


Also from Equation 38, 

F 0 ~ “ ^0 ” “ ^0 ’ 

G 0 = “M 0 G 0 = 0 . 

If we put 


F ' = r F =0 

r 0 r 0 r 0 u ’ 


F " = r 2 F + r n a F 1 - 2 

r o r o r o 0 u o r o - -t q 2 /jl 0 


C o - r o G o - r o - 


G o" = r o 2 G o + r o a o G o = r o 2a o 


(obtained from Equation 39) into Equation 40, we obtain 


F(q) “ 1 - M 0 r 0 2 c 2 q 2 


and 


C(q) - r 0 [c,q + cr g r Q c 2 q 2 ] . 

Differentiating with respect to q and regarding Equation 23, we have 

F'(q) = “ M 0 r 0 2 c jq , 

C'(q) = r 0 [c 0 + %r oCl q] . 

By substituting r 0 q = zr , c 0 = /i g r 2 , and v 0 = o- g r into the above equation, we have: 


F = 1 - c 2 f 0 z 2 , 

G = zt(c, +c 2 Vq z) = r(l-c 3 <f 0 z 3 ) 


( 41 ) 
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because of Equation 31b; and 


F = 7 F' = 


M 0 ZT 
c i A 


G = — G' 


c o + c i ^o z 
S 


A “ c 2 £ 0 z = 


because of Equation 29b. 





> (42) 
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SOLUTION OF THE MAIN EQUATION 

With the formulas of Equations 41 and 42, the problem of calculating the ephemerides 
from given initial values is led back to the solution of the main equation. As was previ- 
ously mentioned, this equation is rational for two cases: (1) the circular orbit, where it 
has the trivial form z = 1; and (2) the parabolic orbit, where it takes the form of a cubic 
equation. In all other cases the main equation is transcendental with the outer form of a 
cubic equation whose coefficients are only slightly dependent upon the unknown z. If the 
intermediate time r is not overly large — and this either is never the case in prac- 
tical applications, or it can be easily avoided — the equation can always be solved by means 
of a rapidly converging iteration process if a suitable approximate solution z = z 0 is avail- 
able. For slightly eccentric orbits, but also for orbits of any given eccentricity, if the in- 
termediate time r is short, success is always possible with z 0 = 1, and for near -parabolic 
ellipses or hyperbolas with z 0 as the solution of the cubic equation (Equation 32). 

For the iteration, the Newton approximation method is preferable; i.e., we try to find 
the zero of the function 


H(z) = z + c 2 7] 0 z 2 + c 3 £ 0 z 3 - 1 . 

For z = z 0 , H ( z 0 ) is a small quantity. Now from Equation 23, 

ZF ( C n + l Z " + 1 ) = C n Z " - 

if we set q = zr/r 0 and indicate that here r represents a constant which is also contained 
in T} 0 and £ 0 . Thus, by differentiation with respect to z, we obtain 


= A . 


clH 

dz 


1 + Cj v 0 z + c 2 £ 0 Z 2 
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Therefore, if z x = z 0 + Sz is an improved value of the unknown, 


H(z 0 ) + Sz • A(z 0 ) + 


The versatile properties of the c-functions also make it possible to transform the 
Taylor series 

0 = H(z 0 ) + Sz • H'(z 0 ) +4 < Sz > 2 * H'^o) + ••• 

(in which the derivatives with respect to z are indicated by primes) into a closed expres- 
sion. If we set 

r o 

z - ~ q - Ai . 

where £ is constant, and if we consider that 


d 3 A , dA 

^T + a2 ^ = 0 


= P n 


follows from Equation 18 and from A = r/r 0 , we can also set 


W = A' --a| = A'/3 


= A'" /3 3 


Therefore 


A'" + (j) 2 A' 
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where (a/j3) 2 = p 0 r 2 = Xo . Furthermore since H' = A, 

H iv + X 0 H" = 0 . 

Now if we develop H(z) around the small approximative value H(z 0 )> = H 0 into a Taylor 
series 

H = H 0 + H 0 ' Sz + -JT H 0 " { Sz ) 2 + -jr H 0 "' ( Sz ) 3 + • • • , 
we can then apply the method used for r(q) and, by introducing the c-functions 

C n = C. [X 0 <5Z) 2 ] 

and eliminating the reciprocal factorials by the recurrence formula (Equation 24), we ob- 
tain the rigid expression: 

H = H 0 + H 0 ' Sz + c 2 H 0 " (Sz) 2 + c 3 H 0 "' (Sz) 3 = 0. (44) 

But, if the index 0 in all cases denotes the fact that the respective quantity for z = z 0 is to 
be taken, 

H o' = A o 

V = A 0 ' 

= c o% + C 1 Vo 

= % ( 1 ~ C 2 Xo Z o) + 2o Z o( 1 “ C 3Xo Z 0 2 ) 

+ Vo “ Xo (H 0 - z o + l) 

H o"' = £o - *o ( H o' - 1) 

^0 - Y o ( A o ~ 

By inserting these values into Equation 44, we can use the approximation of Equation 43 
for the calculation of the small factors (Sz) 2 and (Sz) 3 . This also applies for c 2 andc 3 
in a still greater degree; since the argument x 0 (8z)2 is ver Y small, c 2 = 1/2 andc 3 = 1/6 
could always be used. 
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If the orbit is a parabola, z will be a solution of the cubic equation (Equation 32): 


z + az 2 + bz 3 - 1 , 


where a = 7? 0 /2 and b = £ 0 /6. If we introduce u = z - a/3b in place of the unknown z, 


where 


u 3 + au - /3 , 


3b “ a 2 


27 b 2 - 2a 3 + 9 ab 


Now the only real solution of the above equation is given by the Cardanic formula 


sr * ©* * /i - /(if * (if 


because 3b -a 2 = 1/4 ^2<f 0 - v 0 2 ) ~ t 2 /4 ^2e 0 - a 0 ^j and since the discriminant (P/2) 2 + (a/3) 3 
is always positive, for the parabola (p = 0) , a> = 2m = 2e ; therefore 2 e - a 2 = - a 2 =#= p/r 4 

This quantity is always positive since the orbital parameter is positive (except for the case 
of the straight-line orbit, where p - 0). Therefore a > o. 


If we set 


1 + 27 p 


1 +x ) 2 + |/l ~x 2 l~x ) 
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On the other hand, 


4 a 3 _ 4/3 /a\ 3 _ 1 - x 2 

27 ^ “ W\fi) - X 2 



and therefore 



which is inserted into Equation 45 to yield 


_/3 3m 
a 1 + m + m 2 


with 


m 


y i +x 


This solution process can also be used for near-parabolic orbits (it is 
whether they are elliptical or hyperbolic). If A = 2c 2 r] Q) B = Gc 3 £ 0 , 
g = 3B( A + B ) -A 3 , then the main equation is solved by using the equations: 


X 



m 


( 


1 “X 

1 +x ; 


2g m 

1 +m + m 2 

A 

z — u + 


immaterial 
f = 2B - A 2 , 


■mm i 
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Here the integration will start with A = rj 0 , B = £ 5 , since 2 c 2 and 6 c 3 only vary slightly 
from 1 ; and c 2 and c 3 are determined anew with the argument Xo z 2 in second approximation. 

CALCULATION OF AN EPHEMERIS FROM INITIAL VALUES 

The ephemerides of the coordinates of location and velocity of a celestial body moving 
in an undisturbed conic section orbit around the sun can be calculated easily from given 
initial values (x 0 , y 0 , z Q ) and (x 0 , y 0 , z 0 ) by the following method which is suitable for 
electronic computers. 


Given Quantities 

The given quantities are: x 0 , y 0 , z 0 ; x 0 , y 0 , z 0 ; t = k(t - t 0 ) , where k = 0.0172021. If 
an ephemeris is to be calculated for equally spaced times t n = t 0 +nco, where a> - table in- 
terval in days, n = 0 , 1 , 2 , •••, then the intermediate times are r = nkw 
(unit 1 /k = 58 d - 13244). 

Invariables 

The invariables are: 


r- 2 — y 2 <(- v ^ "f 7 ^ * 

r o x o y o z 0 * 

f o 2ct o = Vo + Vo + z « J o i 

r o 2 % = *o 2 + y 0 2 + *<? ■ 

From these, r 0 , cr Qy co Q are calculated, then 

1 

= 77 ’ e o = “o “ Mo 2 p o = Mo ‘ e o ■ 

0 

= /V* ’ \ = a o T - &o = e o r2 ’ Xo = Po t2 ■ 

Solution of the Main Equation 

The main equation for the transfer from t 0 to t, with r = k(tj - 1 0 ) is 
H = z + c 2 ri 0 z 2 + c 3 £ 0 z 3 -1 = 0. 
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Beginning with z = z 0 (see the comments on page 24), we calculate: 

*■* = Xo z o ; 

i \ 2 

c 2 -2T"4T + 6T-"-: 

1 \ 2 

c 3 = 3T~sT + yr----; 

C 1 = 1 “ ^ C 3 : 

A 0 = i + C, 7, 0 z 0 + C 2 £ 0 z 0 2 ; 

H 0 = Z 0 + c 2 % z 0 2 + c 3 Z 0 3 - l ; 

Sz = , z, = z 0 + Sz , X 2 = Xo z l ■ ■ • ' 

This computation is repeated until the values of z 0 , z u z 2 , •• • no longer change. 


Coordinates for t = tj 

Let z, x 2 , Cj, c 2 , c 3 , A be the values from the last iteration for which H = 0 is ful- 
filled. We then form: 


1 -F = c 2 f 0 z 2 ; 

G = r (l - c 3 z 3 ) ; 


• _ ^o 2 

F " “ c i ~5r • 


1 -G 


1 -F 

— g- : 


x l “ x o = - (1-F)x 0 + Gx 0 : 


= f x 0 ~ (1 


G)i„ 


24 


The appropriate formulas apply for y t , z 1 \y 1 , z r 


Controls 

Once we have calculated the invariables n t , a lf a> 1 , ••• with the new coordinates 
(xj, yj, z t ; ij, yj, ij), then the integral relations 


P_0 

p \ A 2 

0 



must be fulfilled. 

COMMENTS 

In computing an equidistant ephemeris, we have two possibilities: (1) after the first 
step has yielded the coordinates x a , y 1? Zy ; x v y v z 1? we use these elements as initial 
values and t 1 as the initial time in the second step. The computation then yields 
x 2 j ••• z 2 ; the procedure is continued, always with the same intermediate timer. (2) We 
use the same initial values (x 0 , *** z 0 ) and their derived invariables for a larger selection 
of steps with increasing intermediate times. 

Each method has advantages and disadvantages. In the first, since r is always small, 
the iteration for the solution of the main equation can usually be begun with z = i and it 
converges rapidly. A disadvantage, however, is that the rounding-off errors in the series 

X 1 = x 0 + ( X l _X o)- x 2 = x 0 + K _X o) + ( X 2 _X l)> ••• 

accumulate rapidly. Also, every step requires the computation of the velocities (which 
are usually not used in practice), and the invariables. The second method does not have 
these disadvantages: r becomes larger for each step; yet the main equation can be solved 
easily, even if z = 1 does not suffice as the initial hypothesis. In this case we can always 
begin with the z which appeared as a solution of the main equation in the last step. Thus 
if Zj, z 2 , z 3 , ••• z n are the solutions of the main equations for the transfer from t 0 to 
t , t , t • • • t ^ the determination of z is begun with z = z and that of z . with 

L 1 9 2 7 l 3 n ? n n- l 7 *1 



25 


Therefore the second method is generally preferred since it is time-saving and more 
exact. Nevertheless, from time to time, depending upon the size of the table interval se- 
lected — perhaps after every fifth or tenth or fifteenth step — we will transfer to new initial 
elements so that r and 77 , £ , x do not become too large. As long as pr 2 is of moder- 
ate magnitude, the power series for the c-functions converge fairly rapidly. For program- 
ming in electronic computers, the following expression is recommended: 


c „M 



X 2 

(n + 1 ) (n + 2 ) 


X 2 f 

(n + 3 ) !n + 4 ) l 1 


X 2 

(n + 5) (n + 6 ) ’ “ 



} 


which allows the automatic computation from within, where the number of terms in the 
brackets depends upon the magnitude of X 2 . For large values of X 2 this method converges 
too slowly; and the trigonometric form, which is not well suited for electronic computa- 
tion, can be used: 


c 


1 


sin X 

X“ 1 c 



X - s in X 
X 3 


This unsuitability is the principal reason for not allowing the intermediate time to in- 
crease too rapidly. 

The major areas for application of the methods described above are: (1) the improve- 
ment of an orbit which originally was determined using the Laplacian method; and (2) the 
computation of special perturbations. In the former application, the first orbit determina- 
tion always yields the local elements for a given time t 0 . The improvement of these ele- 
ments by the use of actual observations, which can be spread out over a long period of 
time, necessitates the calculation of the coordinates at these times for correlation with 
the observations in question. 

To apply the method to the special perturbations of the right-angle coordinates of a 
planet or of a comet, ephemerides of the undisturbed orbit of the celestial body are re- 
quired, which can be readily calculated by means of the above method. As far as the co- 
ordinates of the perturbing planets are concerned, these values were taken from the year- 
books for the times of the calculation, which ' cald be unsuitable for electronic computation. 
In this case it would be better to use the abcr'e method, by computing (beginning with loca- 
tion and velocity coordinates of the perturbing planet at the initial epoch) an undisturbed 
ephemeris (osculating orbit) of the planet which is sufficiently exact for determining the 
perturbations. The latter method has been utilized in Germany with good results. 
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